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A finite-difference, three-dimensional, incompressible Navier-Stokes formulation 
for calculating the flow through turbopump components is presented. The solution 
method is based on the pseudocompressibility approach and uses an implicit-upwind 
differencing scheme together with the Gauss-Seidel line-relaxation method. Both steady 
and unsteady flow calculations can be performed using the presented algorithm. In 
this paper, the equations are solved in steadily rotating reference frames by using the 
steady-state formulation in order to simulate the flow through a turbopump inducer. 
Eddy viscosity is computed by using the Baldwin- Lomax model. Numerical results are 
compared with experimental measurements and good agreement is found between the 
two. Time-accurate calculations will be reported in future publications. 


INTRODUCTION 

With the advent of supercomputer hardware, as well as fast numerical methods, 
computational fluid dynamics (CFD) has become an essential part of aerospace resear ch 
and design. Numerical studies in incompressible flows show good progress in parallel 
with computational studies in compressible flows. For example, the incompressible flow 
solver developed by Kwak et al. (ref. 1) was extensively used for simul at ing the flow 
through the main engine power-head components of the space shuttle. The redesign 
of the ma i n engine hot-gas manifold of the space shuttle, guided by the computations 
of Chang et al. (ref. 2), illustrates the usefulness of CFD in aerospace research. Since 
the incompressible Navier-Stokes formulation does not yield the pressure field explicitly 
from the equation of state or through the continuity equation, numerical solution of the 
equations requires special attention in order to satisfy the divergence-free constraint on 
the velocity field. 

The most widely used methods that use primitive variables are the fractional- 
step and pseudocompressibility techniques. In the fractional-step method, the auxiliary 
velocity field is solved by using the momentum equations. A Poisson equation for pres- 
sure is then formed by taking the divergence of the momentum equations and by using 
a divergence-free velocity-field constraint. Solving the Poisson equation for pressure 
efficiently in three-dimensional curvilinear coordinates is the most important feature 
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of the fractional-step method (ref. 3). One way to avoid the numerical difficulty im- 
posed by the elliptic nature of the problem is to use the pseudocompressibility method. 
With the pseudocompressibility method, the elliptic-parabolic-type equations are trans- 
formed into hyperbolic-parabolic-type equations. Well-established solution algorithms 
developed for compressible flows can be utilized to solve the resulting equations. 

Steger and Kutler (ref. 4) employed an alternating-direction implicit scheme 
in Chorin’s (ref. 5) pseudocompressibility method. This formulation was extended to 
three-dimensional generalized coordinates by Kwak et al. (ref. 1). Recently, a three- 
dimensional incompressible Navier-Stokes solver (INS3D-LU-SGS) using a lower-upper- 
symmetric Gauss-Seidel algorithm was developed by Yoon and Kwak (ref. 6). This al- 
gorithm is used to calculate the inducer flow of the main-engine turbopump of the space 
shuttle in order to demonstrate the performance of the numerical method (ref. 7). An- 
other effort is reported in reference 8 in which upwind differencing and the Gauss-Seidel 
line-relaxation scheme were used in order to have a robust and fast-converging scheme 
(INS3D-UP). A time-accurate formulation of this algorithm is implemented for incom- 
pressible flows through artificial heart devices with moving boundaries (refe. 8-10). In 
the present study, the steady-state formulation is used in steadily rotating reference 
frames in order to develop a CFD procedure for simulating the flow through the tur- 
bopump components of a liquid rocket engine. 

In the first section, the governing equations and numerical method are described. 
Following that is a presentation of the computed results obtained from the current 
approach. 

The authors would like to thank Mr. Liang-Pang Chang of MCAT Institute for 
providing grid-generation and graphical support for this study. This work is partially 
supported by the NASA Marshall Space Flight Center. Computer time was provided 
by the Numerical Aerodynamic Simulation (NAS) facility at Ames Research Center. 


GOVERNING EQUATIONS AND METHOD OF SOLUTION 


The algorithm used in both steady and unsteady formulations is based on the 
method of pseudocompressibility, which produces a hyperbolic system of equations by 
introducing a time-derivative pressure term into the continuity equation. The result- 
ing incompressible Navier-Stokes equations can be written in a generalized curvilinear 
coordinate system (f, rj, C) as follows: 

§+!(*-*.)+£(*-jw+£<g-g.)-* w 

where Q, and the convective flux vectors E, F, G, are 

V 

u 
v 

w. 

2 




* 


r 0U i 

E = i ZxP + uU £ t u 
J kyP + VU + £ t V 
-tzP + wU + £ t w. 



PV I 

VxP + uV + 7) t U 
VyP + vV + r} t v 
.TJzP + wV + T) t W. 


pw 

Q = }_ CxP + uW + £ t U 
J CyP + vW + ( t V 
-C zP + wW + Ctty. 


Here J, P, p, u, v, and w denote the Jacobian of transformation, the pseudocom- 
pressibility coefficient, pressure, and Cartesian velocity components, respectively. The 
contravariant velocity components U, V, and W are defined as 


U = £xU + £ y v + 

V = ij x u + rjyV + T) z w 
W = C x u + C y v + ( z w 


For an orthogonal-grid assumption, the viscous flux vectors E v , F v , and G v are given 

by 


Fv= Re 1 
Gv = ReJ 
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where Re is the Reynolds number. For simplicity, the viscous terms are written for 
the orthogonal grid and the constant viscosity case. However, the algorithm includes 
the terms necessary for a nonorthogonal-gnd case and for a nonconstant-" viscosity case. 
When the equations are solved in steadily rotating reference frames, the centrifugal 
force and the Coriolis force axe added to the equations of motion as source ter ms In 
this case, the equations are written for relative velocity components in a moving relative 
frame of reference. The source term S is given by 


0 1 

5= 0 
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where ft is the rotational speed. Relative velocity components are written in terms of 
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absolute velocity components u a , v a , and w a as 

u = u a 
v = v a + Qz 
w = w a — Sly 

In other cases, the source term S is set to zero. In the steady-state formulation, the 
time-derivatives are differenced by using the Euler backward formula. The equations 
are solved iteratively in pseudotime until the the solution converges to a steady state. 
In the time- accurate formulation, the time-derivatives in the momentum equations are 
differenced by using a second-order, three-point, backward-difference formula. The 
equations are iterated to convergence in pseudotime for each physical time-step until 
a divergence-free velocity field is obtained. Central differencing is used to numerically 
compute the viscous flux derivatives, because they have well-behaved diffusive terms. 
However, higher-order upwind differencing is employed to compute the convective flux 
derivatives. The reason for using an upwind scheme is to have a more reliable, robust, 
and general-purpose solution algorithm. Chakravarthy and Osher outline the class of 
high-accuracy flux differencing schemes for compressible flow equations (ref. 11). Fol- 
lowing the third-order scheme of reference 11, a fifth-order-accurate, upwind-biased 
stencil is derived by Rai (ref. 12). The upwind differencing used here is the implemen- 
tation of those efforts on incompressible Navier-Stokes equations. The ^-derivative of 
the convective flux E can be written as 

dE ^ [i?i+ 1/2 — £*- 1 / 2 ] 
dt ~ A£ 

The numerical flux E i+ 1 / 2 is defined as follows: 

Ei+ 1/2 = g [^(Q*+i) + E(Qi ) — &+ 1 / 2 ] (2) 

where the fc+ 1/2 is a dissipation term. The order of the scheme is determined by the 
definition of the dissipation term <t>i+i/ 2 - For fa+ 1/2 = 0, the differencing is reduced to 
a second-order central-difference scheme. A first-order upwind flux is defined by 

&+ 1/2 — (A-^i+i /2 — A -^+ 1 / 2 ) (3) 

and a third-order upwind flux is given by 

<f>*+l/2 = ~ 2 ~ ^^i+1/2 + AE i+ 1/2 — AE i+ y^j (4) 

Finally, the following dissipation term creates a fifth-order upwind-biased scheme that 
requires seven stencil points: 

&+ 1/2 = ~ 2 Q ( — 2AJ5£. 3 /2 + 11 A — 6A E? +1 / 2 — 3AEt +3 ^ 2 + 2A E i+5 ^ 2 
- 11A£?- 3/2 + 6A E~ +1/2 + 3A£- 1/2 ) 

where A E ± is the flux difference across positive or negative traveling waves, and is 
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computed as 


A-^h-i/2 = (Q)AQi + i/2 

Here, A ± is the plus (minus) Jacobian matrix. The A operator, and Q are given by 

&Qi+ 1/2 = Qi + 1 — Qi . 

Q — 2^* +1 ^*) 

A special treatment is neccesary next to the boundary when third- and fifth-order 
differencing schemes are used. The differencing near the boundary is 

8E _ [£*=5/2 - £’ 1 = 3 / 2 ] 

dt~ At 

For the flux terms, both at i = 3/2 and at i = 5/2, the first-order terms are used. 
The differencing between the flux at i = 7/2 and the flux at i = 5/2 is increased to 
third-order differencing by using third-order flux terms at both points. This requires 
storing both first- and third-order flux terms at i = 5/2. 

An implicit, delta-law-form approximation to equation (1) after linearization in 
time and the use of approximate Jacobians of the flux differences results in a seven-block 
diagonal matrix equation written as 

B6Qi-ij,k + ASQij'k + C6Q i+ i dtk + D6Qij- lile + ESQ itj+ltk + F6Q itj) k- 1 

+ G6Q i<j)k+1 = R.H.S. {6) 

where 6Q = Q n+1 — Q n and A, B, C, D, E, F, and G are 4x4 block diagonals. 

The solution method of this matrix equation greatly affects the efficiency and 
the robustness of the whole algorithm. An approximate factorization (AF) scheme 
introduces a time-step limitation owing to factorization error (ref. 13). Since the flux- 
difference splitting requires more computer time than central differencing, using an 
unbounded tune-step is important in order to reach the efficiency of a Newton pro- 
cedure. The Gauss-Seidel line-relaxation scheme, which was succesfully em ployed by 
MacCormack (ref. 14), allows an unbounded time-step to be taken. 

In equation (6), the right-hand-side term is computed and stored for the entire 
domain. The Gauss-Seidel line-relaxation procedure is composed of three stages, each 
of which involves a block-tridiagonal inversion in one direction. In the first stage, 6Q is 
solved line-by-line in one direction. Before the block-tridiagonal equation is solved, off- 
tridiagonal terms are multiplied by the current value of 6Q and are then shifted over to 
the right-hand side of the equation. In other words, equation (6) is solved by performing 
a block-tridiagonal inversion in the ^-direction, and Gauss-Seidel sweeps in the rj- and 
(-directions. After the first sweep is completed for the entire domain, a second sweep 
is started in the opposite direction. One forward sweep and one backward sweep Ruffle 
for most of the problems, but the number of sweeps can be increased. The second stage 
is to solve the block-tridiagonal terms in the redirection and to perform backward and 
forward sweeps in the (- and (-directions. The same procedure is repeated in the 
third stage by inverting the block-tridiagonal matrix in the (-direction and treating 
the off-line terms for the (- and redirections in Gauss-Seidel fashion. 
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Implicit boundary conditions are used at all of the boundaries except at the zonal 
and overlapped interface boundaries; zonal and grid-embedding interface boundaries are 
explicitly updated. No-slip boundary conditions are imposed at the solid stationary 
wall where the velocities are zero. The pressure boundary condition is specified such 
that the pressure gradient normal to the wall is zero. 

At the inflow and outflow boundaries, characteristic boundary conditions are 
employed implicitly. It is assumed that the effect of viscous terms at the boundaries 
is negligible. Also, the characteristic equations are approximated in one-dimensional 
space. The primitive variables that are needed at the boundaries are the pressure p and 
the velocity components u, v, and w. The number of positive and negative eigenvalues 
of the Jacobian matrix of the convective flux determines how many variables should 
be specified at the boundaries. If the flow is in the positive ^-direction then three of 
the eigenvalues are positive and one of them is negative. Therefore, there will be three 
characteristic waves traveling downstream and one characteric wave traveling upstr eam 
In that case, the exit boundary receives information about the three variables by means 
of the characteristics traveling from the interior of the domain. Hence, only one variable 
should be specified at the outflow boundary. However, the inflow boundary is informed 
by only one characteristic traveling from the interior region. Therefore, three variables 
should be specified at the inflow boundary. Details of the numerical method are given 
in reference 8. 


COMPUTED RESULTS 

The flow field through a turbopump inducer is solved as a benchmark prob- 
lem in order to validate the CFD procedure for turbomachinery applications. In this 
section, results obtained for the Rocketdyne inducer shown in figure 1 are presented. 
The inducer geometry was developed and experimentally studied by the Rocketdyne 
Division of Rockwell International. The design flow is 2,236 gal / min with a design 
speed of 3,600 rpm. In the computational study, tip-leakage effects are included with a 
tip clearance of 0.008 inch. The problem was nondimensionalized with a tip diameter 
of 6.0 inches and an average inflow velocity of 28.3 ft /sec. The Reynolds number for 
this calculation was 191,800 per inch. The upstream section of the inducer was taken 
as a 2-tip-diameter-long straight channel, as shown in figure 1. The bull-nose of the 
inducer was treated as a rotating wall, and the cavity section was neglected. However, 
this region can be included by using an additional zone. An H-H grid topology with 
d imen s i ons of 187 x 27 x 35 was used. A partial view of the surface grid is shown 
in figure 2. An H-type surface grid was generated for each surface using an elliptic 
grid generator. The interior region of the three-dimensional grid was filled using an 
algebraic solver coupled with an elliptic smoother. In the straight channel, the grid 
was generated for one-sixth of the cross section of the tube. This grid was extended to 
the outflow section of the inducer between the blades. Periodic boundary conditions 
were used at the endpoints in the rotational direction. 

At the inflow and outflow boundaries, characteristic boundary conditions were 
employed. At the inflow, v and w velocity components were specified as zero, and 
the total pressure was specified as constant. Axial velocity and static pressure were 
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calculated from the characteristic relation and the total pressure relation. At the 
outflow, static pressure was specified and the velocity components were computed from 
the characteristics propagating from the interior region. The flow was taken at rest 
initially and the inducer was fully rotated impulsively. The solution was considered 
converged when the maximum residual dropped at least four orders of mag nitude 
This was obtained in less than 500 iterations. The computer time required per grid 
point per iteration was about 1.4 x 10 -4 sec. 

Figure 3 illustrates the planes where the experimental measurements were taken 
by Rocketdyne. Axial and tangential velocity components and the flow angle were 
measured in planes A, B, C, and D at various circular arcs from the hub to the tip region. 
At each plane, the comparison between experimental measurements and numerical 
results along three of the circular arcs is presented in this paper. A total velocity and 
a flow angle are compared with the experimental data. The total velocity has only a 
tangential and an axial velocity component. The radial velocity component was not 
measured in the experiment. 

Figures 4 through 7 show relative total velocities and relative flow angles as a 
function of circumferential angle in planes A, B, C, and D, respectively. The circum- 
ferential angle increases from the suction side to the pressure side. The comparison of 
computed and experimental results is generally good all the way from the hub to tip re- 
gion. The difference between the experimental data and the numerical results is about 
5%-8% in velocity. In all planes, the hub and tip regions show the biggest discrepan- 
cies. This may be a result of the relatively coarse grid used for the boundary layer. 
In the computational study, the Baldwin-Lomax algebraic turbulence model is used to 
determine the eddy viscosity. The comparison shows that the solution algorithm does 
a good job with an algebraic turbulence model. Implementation of the one-equation 
model of Baldwin and Barth (ref. 15) is under way for the present algorithm. 

The reason for higher-order turbulence modeling is the comparisons obtained in 
plane D, in which the wake region is not predicted accurately (fig. 7). Another advan- 
tage of the one-equation model is that there is no need to define a length-scale explicitly. 
Near the tip-clearance region, the difference between experimental measurements and 
numerical results is noticeably larger than the error in other regions. This is due to 
the lack of grid resolution in the tip-clearance region. In the grid-refinement study, the 
number of grid points in the tip-clearance region was increased from four to nine. In the 
coarse-grid computation there is one overlapped grid point in the rotational direction 
to ensure periodic boundary conditions. In the fine grid, three additional zones were 
added in radial direction. The results with the one-equation model and the results from 
the grid-refinement study will be published later. 

Figure 8 shows the surface of the inducer color-coded by nondimensionalized 
pressure. The pressure gradients across the blades caused by the action of centrifugal 
force and the pressure rise from inflow to outflow are illustrated. This pressure rise along 
the inducer can also be seen in figure 9. Velocity vectors are plotted in the meridional 
plane and the vectors are colored by the static pressure. The existing solution proce- 
dure can be applied to the same configuration under off-design conditions. The massive 
separation, which may block the fuel supply, can be detected in the numerical study. 
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This separation will be the subject of our future work which is intended to produce a 
pre-design and post-design engineering tool in challenging turbo machin ery applications. 

CONCLUDING REMARKS 

An efficient and robust solution procedure was developed and validated for three- 
dimensional turbopump applications. Numerical simulations of the flow through a 
Rocketdyne inducer were successfully carried out by using CFD te chni ques for solv- 
ing viscous incompressible Navier-Stokes equations with the source terms in steadily 
rotating reference frames. The method of pseudocompressibility with a higher-order 
accurate upwind differencing and the Gauss-Seidel line-relaxation scheme provide fast 
convergence and robustness. Results in the form of relative total velocity and relative 
flow angle in four planes are presented. Numerical results compare fairly well with 
experimental data. 


8 



REFERENCES 


1. Kwak, D.; Chang, J. L. C.; Shanks, S. P.; and Chakravarthy, S.: A Three- 

Dimensional Incompressible Navier-Stokes Flow Solver Using Primitive Vari- 
ables. AIAA J., vol. 24, no. 3, 1977, pp. 390-396. 

2. Chang, J. L. C.; Kwak, D.; Rogers, S. E.; and Yang, R.-J.: Numerical Simulation 

Methods of Incompressible Flows and an Application to the Space Shuttle 
Main Engine. Int. J. Num. Meth. Fluids, vol. 8, 1988, pp. 1241-1268. 

3. Rosenfeld, M.; Kwak, D.; and Vinokur, M.: A Fractional Step Solution Method 

for the Unsteady Incompressible Navier-Stokes Equations in Generalized Co- 
ordinate Systems. Journal of Computational Physics, 1991, in press. 

4. Steger, J. L.; and Kutler, P.: Implicit Finite-Difference Procedures for the Com- 

putation of Vortex Wakes. AIAA J., vol. 15, no. 4, Apr. 1977, pp. 581-590. 

5. Chorin, A. J.: A Numerical Method for Solving Incompressible Viscous Flow 

Problems. J. Comput. Phys., vol. 2, 1967, pp. 12-26. 

6. Yoon, S.; and Kwak, D.: Three-Dimensional Incompressible Navier-Stokes Solver 

Using Lower-Upper Symmetric-Gauss-Seidel Algorithm. AIAA J., vol. 29, 
no. 4, 1991, pp. 874-875. 

7. Yoon, S.; and Kwak, D.: Implict Methods for the Navier-Stokes Equations. Com- 

put. Sys. Eng., vol. 1, nos. 2-4, 1990, pp. 535-547. 

8. Rogers, S. E.; Kwak, D.; and Kiris, C.: Numerical Solution of the Incompressible 

Navier-Stokes Equations for Steady and Time-Dependent Problems. AIAA J., 
vol. 29, no. 4, 1991, pp. 603-610. 

9. Kins, C.; Chang, I.; Rogers, S. E.; and Kwak, D.: Numerical Simulation of 

the the Incompressible Internal Flow through a Tilt ing Disk Valve. AIAA 
Paper 90-0682, 1990. 

10. Kiris, C.; Rogers, S. E.; Kwak, D.; and Chang, I.: Computation of Incompress- 

ible Viscous Flows with Moving Boundaries. Proceedings of the International 
Symposium on Biofluiddynamics, Seattle, Wash., July 1991. 

11. Chakravarthy, S. R.; and Osher, S.: A New Class of High Accuracy TVD Schemes 

for Hyperbolic Conservation Laws. AIAA Paper 85-0363, 1985. 

12. Raj, M. M.: Unsteady Three-Dimensional Navier-Stokes Simulations of Turbine 

Rotor-Stator Interaction. AIAA Paper 87-2058, 1987. 

13. Beam, R. M.; and Warming, R F.: An Implicit Finite-Difference Algorithm for 

Hyperbolic Systems in Conservation-Law Form. J. Comput. Phys., vol. 22, 
1976, pp. 87-109. 

14. MacCormack, R. W.: Current Status of Numerical Solutions of the Navier-Stokes 

Equations. AIAA Paper 85-0032, 1985. 

15. Baldwin, B. S.; and Barth, T. J.: A One-Equation Turbulence Transport Model 

for High Reynolds Number Wall-Bounded Flows. AIAA Paper 91-0610, 1991. 


9 





Figure 2. Surface grid for Rocketdyne turbopump inducer. 
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Figure 3. Schematic representation of planes where experimental data are available. 
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Figure 5. Comparison of relative total velocity and relative flow angle in plane B. 
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Figure 6. Comparison of relative total velocity and relative flow angle in C. 
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Figure 8. Surface pressure for Rocketdyne inducer. 
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Figure 9. Velocity vectors colored by pressure on the meridional plane of the inducer. 
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